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The stability of the dynamical states characterized by a uniform firing rate {splay states) is 
' analyzed in a network of globally coupled leaky integrate-and-fire neurons. This is done by reducing 

' the set of difi'erential equations to a map that is investigated in the limit of large network size. We 

, show that the stability of the splay state depends crucially on the ratio between the pulse-width 

' and the inter-spike interval. More precisely, the spectrum of Floquet exponents turns out to consist 

of three components: (i) one that coincides with the predictions of the mean-field analysis [Abbott- 
van Vreesvijk, 1993] ; (ii) a component measuring the instability of "finite-frequency" modes; (iii) a 
number of "isolated" eigenvalues that are connected to the characteristics of the single pulse and 
may give rise to strong instabilities (the Floquet exponent being proportional to the network size). 
Finally, as a side result, we find that the splay state can be stable even for inhibitory coupling. 
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I. INTRODUCTION 



, Understanding the mechanisms of information processing in the brain can be tackled by analyzing the dynamical 
properties of neural network models. Nonetheless, approaching the problem in full generality is a formidable task, 
since it requires taking into account, (i) the role of the topology of the connections, (ii) the dynamics of the connection 
'"^ > themselves as a representation of the synaptic plasticity, (iii) the internal dynamics of the model neurons, which can 
^ ' depend on the number of ionic channels and other variables and parameters, (iv) the diversity among neurons and 
O , connections, (v) the unavoidable presence of noise. On the other hand, we can conjecture that at least some basic 
mechanisms are quite robust and may depend on a few ingredients, only. In fact, even simple models made of globally 
coupled identical units exhibit interesting dynamical properties, far from being completely interpreted. For instance, 
the stability of their steady states is still a debated problem. In fact, one can find claims that the splay state [l|, 
characterized by a uniform spiking rate of the neurons, is stable only in the presence of excitatory coupling 0, and 
yet stable splay states have been found also in networks with fully inhibitory coupling [5]. The standard approach 
for determining the stability properties of such simple models is based on the mean field approximation. This allows 
to obtain the spectrum of eigenvalues associated with the stability matrix in the thermodynamic limit N ^ 00, N 
denoting the number of neurons. It is far from obvious if and how such results can be extended to the stability 
problem of large but finite networks. This is all more important if we consider that in the thermodynamic limit the 
^ ^ splay state has been shown to be marginally stable for finite pulse-width, while it has been found to be strictly stable 
[ for (5-like pulses [3, Moreover, we have discovered that in the limit of vanishing pulse-width, the formula derived 
in does not coincide with the result of direct simulations. 

We want to point out that assessing the linear stability of the splay state is just a preliminary step towards a 
complete undertanding of the dynamical properties of neural networks. Nonetheless, even the accomplishment of 
' this "simple" task is going to reveal unexpected subtleties, which should be taken into account for pursuing further 
progresses. In particular, in this paper we introduce a specific formalism for determining the stability of the splay 
state in (in)finite ensembles of neurons. The approach is not entirely new, as it is based on the introduction of a 
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Poincare section, which transforms the original dynamical system into a map connecting dynamical configurations of 
the neural network after consecutive neural pulses (spikes) . A similar technique was used in [s*] for obtaining analytic 
upper bounds of the maximum stability exponent of disordered neural networks with global inhibition. In 6] such a 
technique was exploited to setup a general approach for the investigation of various dynamical regimes, including the 
splay state. The analysis was performed by considering the dynamical transformation describing consecutive spikes 
of the same neuron: due to its complicacy such a transformation does not allow to proceed much beyond formal 
statements. In this paper we have taken advantage from the idea of constructing a map between consecutive spikes 
of whatever neurons, combined with a suibtale shift of the neuron labels. In this way the computational complexity 
is significantly reduced with respect to 0] and we are able to obtain analytic expressions also for large finite values of 
N. 

One important consequence of our analysis is that one must be extremely careful when dealing with the stability 
problem genereated by finite pulse-width in large networks. Actually, even if the leading correction to the mean-field 
approximation is found to be 0{1/N), we show that higher order corrections, up to 0{1/N^) , have to be taken into 
account for determining the correct stability properties, even for arbitrarily large values of N. 

Another important result of our study is that the thermodynamic limit does not commute with the zero pulse-width 
limit. In order to clarify this issue we consider the stability problem in the presence of pulses, whose width is inversely 
proportional to N . Although this recipe may appear a mathematical trick, it allows to point out a crucial aspect of 
the problem at hand: when both N and the spiking rate are large, their ratio is the only relevant stability parameter. 
This implies also that the stability of a network subject to fixed small pulse-widths may qualitatively change when 
N is increased, precisely because the above mentioned ratio varies. 

In this paper we also clarify the basic question whether a given network of pulse-coupled neurons exhibits a finite 
stability or if it is "marginally" stable. A general formula obtained in [3| indicates that the amplitude of the eigenvalues 
of the stability spectrum of the splay state converges to zero from below as This implies that infinitely large 

networks exhibit an arbitrarily large number of eigenvalues arbitrarily close to zero. Accordingly, one should conclude 
that marginal stability is a generic feature of infinite networks. This seems to contrast with the result of Ref. [sj, 
where numerical evidence of stability in the thermodynamic limit was found for 5-like pulses. Our analysis clarifies 
that the Abbott- Van Vreeswijk's formula [3] applies to the case of finite pulse-width, while in the case of vanishing 
pulse-width, stability is maintained also in the thermodynamic limit. The following argument provides a heuristic 
explanation of the mechanism generating such different stability conditions. In networks of pulse-coupled oscillators 
the same pulse acts differently on the various oscillators, because they can be located in different regions of their phase 
space, when the pulse is received. If the states of two neurons are close to each other in phase space (e.g., the neurons 
exhibit almost equal values of their membrane potentials), they will experience only slightly different effects from the 
same pulse of finite width. This amounts to say that the two neurons are very weakly coupled, thus yielding very small 
values of the Floquet stability exponent. The same is no longer true when very small (infinitesimal) pulse-widths are 
considered, since the coupling field strongly oscillates also over "microscopic" time scales, even for large values of N . 
In other words, the "roughness" of the coupling is able to strongly lock the neurons, thereby yielding a finite stability 
of the splay state. 

By this argument one can also understand the limitations of the mean- field approach proposed in 0. In fact, 
once the neurons are ordered according to their instantaneous potential, we have verified that in many cases the 
stability of the state is controlled by the stability of "high-frequency" modes. A typical example of such a mode is the 
perturbations associated with the forward (backward) shift of the even (odd) neurons (see Section IV). These modes 
are by definition neglected in the mean-field approach _2:]. 

In Sec. II we describe the model of leaky integrate-and-fire neurons considered in this paper and we reduce it in 
full generality to an event-driven map. Moreover, also the formalism for the linear stability analysis is introduced. In 
Sec. Ill, we discuss the stability of the splay state in the case of finite pulse-width. Sec. IV is devoted to the analsysis 
of vanishing (l/iV) pulse- widths. A brief summary of the results and future perspectives are reported in Sec. V. 

II. THE MODEL 

We consider a simple model of a network of identical leaky integrate-and-fire neurons, whose interaction is mediated 
by a field E{t) . This field can be viewed as the result of the linear superposition of single-neuron pulses Es{t), emitted 
when the membrane potential reaches a given threshold value. In agreement with Ref. 01, we assume that the time 
dependence of a single pulse is given by the relation Es{t) — a^te~°'* , where 1/a is the pulse-width. Accordingly, in 
between consecutive spikes, the field evolution is described by the differential equation. 



E(t) + 2aE{t) + a^E{t) = 



(1) 



Moreover, the effect of the pulse emitted at time to amounts to a discontinuous change of the derivative, 



The dynamics of the membrane potential Xi{t) of the i-th neuron is determined by the differential equation 

Xi = a ~ rjXi + gE{t) , e (-00, 1) , (2) 

where the parameter g gauges the strength of the coupling with the pulse field E{t). The membrane potential evolves 
according to ^ until it reaches its threshold value Xi = 1 , when it is instantaneously reset to the value Xi = . Such 
a condition amounts to representing the discharge mechanism in neurons as a strong nonlinear effect associated with a 
discontinuity in the membrane potential. One should point out that this model is written from the very beginning in 
adimensional rescaled units (for a comparison with physical scales see the discussion in ref. Q ) . Without prejudice 
of generality one can further simplify the equations by fixing the scale of time in units of the parameter 77"^ = 1. 
Consistently, one should formally rescale also the parameter a to a/77. 

A. Event-driven map 

By integrating Eq. ([T]) between successive pulses, one can reduce the continous-time evolution of the interaction 
field E{t) to the discrete map, 

E{n + 1) ^ E{n)e-°'^ + NQ{n)Te~°'^ (3a) 
g(n + l) = g(n)e-"- + ^ , (3b) 

where r is the interspike time interval and we have introduced the new variable Q := {aE + E)/N . 

Accordingly, also the differential equations ([2|) can be integrated by taking into account this map- like representation 
of the field evolution and reduced to an event-driven map for the membrane potential Xi , 

Xi{n + 1) = x.,in)e^^ + a (l - e"^) + gF{n) i = l,...,N . (4) 

where F is 

n„).£::^(^,„, + ^^)-f£Z_„<3,„,. 

a — 1 \ a — 1/ [a — 1) 

If the integer m labels of the closest-to-threshold neuron, the interspike time interval r is obtained by imposing 
the condition, x„i{n + 1) = 1 , 



T(n) = In 



Xmjn) - a 
1 - gF{n) - 



(6) 



Since in networks of identical neurons the order of the potentials Xi is preserved, it is convenient to introduce a 
comoving spatial index frame, j{n + 1) = j{n) — 1 , where the next firing neuron is always labelled by j(n) = 1 and 
reset to j(n + 1) — N after the firing event. Accordingly, the map for the membrane potentials transforms into 

Xj-i{n + I) ^ Xj{n)e^^ + 1 - xi{n)e^'' j = 1, . . . , TV - 1 , (7) 

with the boundary condition xn — 0. The set of Eqs. (|3l6l7p defines a discrete-time mapping (notice that in the 
comoving frame the label m in Eq. © is always set equal to 1) , that is fully equivalent to the original set of ordinary 
differential equations. Altogether, it turns out that a network of N identical neurons is described by -|- 1 equations: 
Two of them account for the dynamics of E(n), while the remaining A^ — 1 equations describe the evolution of the 
neurons {xn is no longer a variable being, by definition, equal to 0). Accordingly we see that, at variance with 
Ref. Q where no field dynamics was explicitely introduced, the model has a finite dimension. 

In this framework, the periodic splay state reduces to a fixed point that satisfies the following conditions, 

T{n) = ^ , (8a) 

E{n) EE E , Q{n) ee Q , (8b) 
= x,e-^^^ + 1 - iie-^/^ , (8c) 



where T is the time elapsed between two consecutive spike emissions of the same neuron. A simple calculation yields, 
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E^TQ 
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The solution of Eq. ([8t) involves a geometric series that, together with the boundary condition ijv = , leads to a 
transcendental equation for the period T . Not surprisingly, the result turns out to be independent of the pulse-width 
a . For simplicity, we report only the leading terms for ^ 1 , 

aT 



T 



T = In 



9 



(9b) 



.(a-l)T + 5j 

By assuming that the uncoupled neurons are in the repetitive firing regime, i.e. by taking a > 1, the period T is 
well defined in the excitatory case (g > 0) only for g < 1 (T — > 0, when g approaches 1), while in the inhibitory case 
{g < 0), a meaningful solution exists for any coupling strength (T ^ oo for g —>■ —oo). 



B. Linear stability 

The main complication of the linear stability analysis stems from the implicit dependence of r on the relevant 
variables (see Eq. ©, where m = 1 in the comoving frame we have adopted). By linearizing Eq. ([H]), we obtain 

5r(ri) = TxSxi{n) + teSEIu) + rgSQin) ; 

where Tx := dr/dxi and analogous definitions hold for te and tq. The linear stability of the splay state is thus 
determined by the dynamics of the perturbations dxj, 5E, and 5Q. This is expressed by the following set of linear 
equations 

5xj^i{n + 1) = e-'^'^[5xj{n) ~ bxx(n)\ + e~^/^(ii - £j)(5r(n) , (10a) 
bE{n + 1) = e-'^^l^hE(n) + Te-"^/^(5Q(n) - (aE - A^Qe""^/^) bT{n) , (10b) 

bQ(n + 1) = e-"^/^5Q(n) - aQe--^^ 1^ 6T(n) , (10c) 

which have been linearized around the fixed point ([5]) . The boundary conditon x^q = imposed by the comoving 
frame yields Sxn — . In practice, the stability problem is solved by computing the Floquet spectrum of multipliers 
{/ife} , fc = 1, . . . , A + 1 associated with the eigenvalue problem of the set of linear equations (fTU|) . In general the 
explicit computation has to be performed numerically. 

For the sake of clarity, it is first convenient to discuss the trivial case of vanishing interaction, i.e. g — . After 
simple calculations the eigenvalues turn out to be fik — exp(i(pfc), where (pk = ^jT' ^ ^ 1, . . . , A — 1, and /ijv = 
(iN+i = exp(— aT/A^) . The last two exponents concern the dynamics of the coupling field E{t), whose decay is ruled 
by the time scale a"^ . 

As soon as the coupling is on, small amplitude fluctuations ^ 0{g/N) affect the neuron dynamics and the spectrum 
of Floquet multipliers takes the general form 

^,^e^^'=e^(^^+-^)/^, = A: = l,...,Ar-l, (11) 

where, Xk and ujk are the real and imaginary parts of the Floquet exponents. As an example, in Fig. [T] we show the 
spectrum of the Floquet multipliers of the splay state for excitatory coupling {g > 0) and finite values of A^ . The 
multipliers with k — 1, . . . , N ~ 1 are very close to the unit circle, while the two isolated multipliers fj,N and fiN+i lay 
very close to the real axis inside the unit circle. We want to point out that already for g/N w 0(10^^) the multipliers 
of the coupled case can be viewed as small "perturbations" of the uncoupled one. This observation supports the 
perturbative approach that we are going to discuss in the following sections. 

Before accomplishing this task, it remains to comment that the variable ipk plays the same role of the wavenumber 
in the linear stability analysis of spatially extended systems, so that we can say that Xk characterizes the stability 
of the fc-th mode. In the following analysis it is convenient to distinguish between modes characterized by ifk ~ 0, 
mod (27r) + 0{1/N), and the other modes. Actually they identify two spectral components, that require a different 
mathematical treatment. The first component corresponds to the condition — 1|| A^~^ and is referred to as 
long wavelengths (LWs), while the second one corresponds to — 1|| ~ 0{1) and is referred to as, short wavelengths 
(SWs). 
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FIG. 1: (Color online) Unitary circle and exact Floquet multipliers spectra for the complete maps for A'' = 20 (red circles) and 
N = 10 (blue crosses) for excitatory coupling. The parameters are a = 3.0, g = 0.4, and a = 30.0. 



III. FINITE PULSE-WIDTH 



In this section we investigate the stabihty problem of the splay state for networks subject to pulses with finite a (i.e. 
independent of the system-size) for large values of N. This is also the setup studied in [2^1 by a mean field analysis. 
By retaining all the terms up to the order the event-driven map (see Eqs.® and ^ ) simphfies to the set of 

+ 1 equations 

E{n + l) = {l-aT)E{n)+NQ{n)T , (12a) 

Q{n+l) = il-aT)Q{n) + — , (12b) 
Xj-i{n + 1) = (1 - T)xj{n) + l-xi{n)+T , (12c) 
where j = 1, . . . , N — I. Here the expression of interspike interval ^ simplifies to 

nn) = . . jp/ ^ ■ 13) 

a — 1 + gtj[n) 

The periodic solution for the pulse-field becomes E = T^^ , and Q — a/NT , while xj and the period T are still given 

by Eq. The Floquet eigenvalue spectrum Hk , k — l,...,iV+l, can be obtained by solving the set of + 1 linear 
equations 

HkSxj^i = (1 - T/N)Sxj - Sxi + (1 - Xj )St , (14a) 

likSE = {l-aT/N)SE + TSQ, (14b) 

fikSQ = {l-aT/N)SQ~—ST. (14c) 

An explicit expression for 5t is obtained by evaluating the derivatives of Eq. (|13p 

St ~ SE Sxi . 

N g 

With this substitution in Eqs. p^ we conclude that the eigenvalue problem amounts to finding the A^-|- 1 roots /i^ of 
the associated polynomial. A partial simplification of the problem can be obtained by extracting from Eas. (|14l b.c) 
the dependence of 5t directly in terms of the perturbation of the potential of the next-to-threshold neuron 5xi^ 

St = KSxi 

where 

a^gT 



K = 



-(a-l + 5/T) + 



N^{^lk - l + aT/Nf 



(15) 



Finally, by substituting this expression into into Eq. (|14b.). we obtain a closed set of equations for the perturbations 
of the membrane potentials, 

Aifefaj_i = (1 - T/N)5x.j + {K - l)5xi - KijSxi . (16) 
After imposing the boundary condition, 5xn = , Eq. (|16p reduces to the following eigenvalue equation, 

^r^e^ ^ K{a + g/T)'—^ - [K{a 1 + g/T) + 1] \ J''^ ' , (17) 

where if is a function of /i^ (see Eq. (fT5|)) . 

In order to solve analytically the above equation, it is necessary to distinguish between long and short wavelengths. 



Let us consider the modes for which ||/ifc — 1| 
simplify the notation we define 



Long wavelengths 

, or, equivalently, ifk ~ 0, mod {2tt) + 0{1/N) . In order to 



At leading order, K writes 



K 



-ia-l + g/T) + 



T(Afc + a)2_ 

Moreover, making use of the relation E — T~^, the eigenvalue equation (|17p can be approximated as 

^)(Afc + a)2(Afc + l)=Afe(e^-e-^''^) , ||Afc||^0. 



a'^g 



(18) 



This equation coincides with that one derived from the mean field analysis 0, except for the missing ||Afc|| = 0, 
which corresponds to the trivial zero Floquet exponent of the continuous time evolution Q and disappears in the 
discrete-time dynamics. We want to stress again that, despite Eq. yields + 1 roots, it provides a suitable 
approximation only for those eigenvalues which satisy the relation — 1|| ^ A^^^ . 



B. Short wavelengths 

The second component of the spectrum (fTT|) is obtained for — 1|| ^ 0{1). In this case K has the simple form 

K = -ia-l+g/T)-' . 

Upon introducing into Eq. I\l7\i the explicit form of the periodic solution ^}p) , the spectrum simplifies considerably, 
namely 

.f = e-^^±^^l, (19) 

i.e., it coincides with a fully degenerate Floquet spectrum, 

ujk = 0, Afe = (20) 

Notice that this approximation holds only for those eigenvalues such that ||/ife — 1|| ~ 0{1), i.e. those representing 
the large maiority of the spectrum, except those laying close to the point (1,0), where the unit circle intercepts the 
real axis (see Fig. [T]) . 

For what concerns the isolated eigenvalues and fiN+i one can easily argue that for finite A^ and finite pulse- 
widths they are always contained inside the unit circle close to the real axis and they approch the point (1,0) from 
the left as A — > oo . Accordingly, they can at most contribute to the marginal stability of the dynamics. 
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FIG. 2: Phase diagram for the stabihty of the splay state in a neural network with excitatory coupling acting through finite 
pulse-width. The solid line separating the stable from the unstable regions in the {g, Q)-plane has been derived from the 
analytic formula of the Floquet spectrum (|18p with a = 1.3. Please notice that in this context stable refers to finite A^, for 
infinite systems the stability will become marginal. 

C. Phase diagram and finite-size corrections 

According to the above results one can conclude that in the hmit N ^ oo the onset of instabihties is determined 
by the Floquet exponents associated with LCs (see Eq. (HH)). It is, therefore, not surprising to discover that our 
result coincides with the predictions of the mean-field analysis reported in [3| . We also confirm that the splay state 
is always unstable for inhibitory coupling {g < 0). In the case of excitatory coupling the mean-field analysis predicts 
stability of the the splay state for a < ac{g,a) , where the ac{g,a) is the critical line separating the stable from the 
unstable region (see Fig. [5])) . By including the role of SWs we can conclude that in the limit N ^ oo the splay state 
can be at most marginally stable for a < ac{g,a) . In particular, ac diverges to +oo for g ^ I (see Fig. [2]). Our 
analysis confirms also that the splay state becomes always unstable via a Hopf-bifurcation, that gives rise to collective 
periodic states 0, [1] . For 5 > 1 no splay state can exist and this can be viewed as a drawback of the model: for too 
strong excitatory coupling, no stationary regime can be sustained, since the evolution steadily accelerates. 

The perfect degeneracy of the zero Floquet exponents associated to SWs sheds doubts on the effective stability 
properties of large but finite networks, since such modes turn out to be marginally stable. Therefore, we have decided 
to solve numerically the eigenvalue equations (dU]) for different system sizes and for values of the parameters g and a, 
which correspond to marginally stable splay states in the limit N 00. The results plotted in Fig.[3|) show that the 
splay state is strictly stable in finite lattices and that the maximum Floquet exponent approaches zero from below as 
This implies that a 1/N perturbation theory, like the one developed in this section, cannot account for such 
deviations and even a second-order approximation scheme cannot capture the instabilities of the original model. This 
is confirmed in Fig. HI where the Floquet spectra obtained from first- and second-order approximations are seen to 
yield an unstable splay state, even though the numerical solution of the stability problem indicates that the finite- A'^ 
model is stable. 

As the event-driven map is the standard approach adopted to simulate this type of networks, an important conse- 
quence of our analysis is that one should be careful enough to consider at least third order approximations in the 1 /N 
perturbative expansion, otherwise the resulting approximate equations may fail to reproduce the correct asymptotic 
dynamics. 

IV. VANISHING PULSE-WIDTH 

In 0] , we have investigated the synchronization properties of a network of leaky integrate-and-fire neurons in the 
case of (5-like pulses. One might expect that the stability properties of such a network can be determined by taking 
the limit of vanishing pulse-width in the formulae obtained in Ref . [3] . However, the lack of agreement with numerics 
suggests that the N 00 limit does not commute with the zero pulse-width limit. In order to clarify this issue, we 
introduce a new setup, where the pulse- width is rescaled with the network size TV as . This amounts to impose 




10" 

l\l 

10 
10" 
10" 
10" 



-3 



: 1 1 1—^ — 1 — 1 1 1 1 1 1 1 — 1 — 1 1 1 1 


1 1 1 — 1 — 1 — 1 1 1. 


- \ 
■ \ 
• ■ \ 


(a): 


■ ■ \ 




\ 


■ 






1 1 1 1 1 ' ' ' 1 1 1 1 1 1 ' ' ' 


\ 



10^ 



10' 



10' 



\x10 



10" 




-600 -400 



-200 
(pkX N 



200 400 600 



FIG. 3: (Color online) (a) Log-log plot of the absolute values of the the Floquet exponents Afc, ordered from the largest to the 
smallest as a function of the index k = 1, N for — 100, 200, 400. The dashed line has slope -2. (b) The Floquet exponent 
as a function of the rescaled phase ifiN, for — 100 (black circles) and A'^ = 200 (red crosses) . In both pictures the parameter 
values are a — 3.0, g = 0.4, and a = 30.0. 
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FIG. 4: (Color online) Floquet exponents Xkif) as a function of the phase ipk for finite pulse-width a = 3.0 and finite neuron 
numbers A'^ = 500 in the case of excitatory coupling g = 0.4. The filled circles represent the exact result for finite A'', while the 
red empty squares and the blue empty triangles refer to approximated results correct up to the first or second order in 1/A'', 
respectively. The parameters here are a = 0.3, rj = 0.1. 



a decay rate proportional to N 

a := (3N. 

It is important to notice that in this context we have to deal with two time scales: (i) a scale of order 0(1) that 
corresponds to the evolution of the the membrane potential xj; (ii) a scale of order ~ that corresponds to 

the field relaxation. 

At leading order, the expression ([5]) for the auxiliary quantity F{n) simplifies to 

F{n) ^ ^ [mn] (1 - e-^^^) + Q{n) (l - e'^f'^ - N^r e"^^-)] . 
The splay state corresponds to the fixed-point 



while the derivatives of the interspike interval ^ now read 



Finally, the eigenvalue spectrum is defined by the set of linear equations 

/ifcfaj_i = {l-TN^^)dxj -dxi+N-\l-S;j)ST, (21a) 
HkSE = e^/^^SE + Te-'^^SQ - N (/3E - Qe~''^^) 6t , (21b) 



= e-l^^dQ - NPQe-l^^dT. (21c) 

We repeat the same formal analysis illustrated in Sec. IIIIl by first solving the eigenvalue equations for the field 
perturbations Sxj . After some tedious calculations we obtain St = Kdxi , where 

K = J . (22) 

(A^fce^^ - 1)2(1 - a - .gii;) + /ifc/32rge/5i^ 



By substituting this expression into Eq. (j2ip and by imposing the boundary condition Sxn — 0, we obtain a closed 
equation for the spectrum {jik} , k = l,...,iV+l, 



where C := {g + aT)/T and K is still a function of ^fc. Proceeding along the previous guidelines, it is necessary to 
separately discuss the behaviour of long and short wavelengths. 



A. Long wavelengths 

As well as in Sec. IIII Al we consider the spectral component characterized by \\^k ^ 1|| ~ ■ In this case K 
assumes the simple expression 

K = -^. 

a - 1 

By neglecting 1/A^ terms, the spectrum satisfies the equation, 

(1 - e^^O (1 + A.T^) = (l - e^(^^+i)) , (24) 

where A^ := NT~^ In/ife . This result agrees again with the mean field analysis [2| in the limit a ^ iV ^ 1 . Let us 
remark that, despite Eq. pi)) yields + 1 eigenvalues, it provides a good approximation only for those such that 

\\^lu-l\\^N~\ 



B. Short wavelengths 



For ||Mfc-l|| 0(l) we can replace the term (/ifee'^^-1) with (e*'^''+^^-l), in Eq. jH]), thus obtaining AT = K{^k) , 
for fc = 1, . . . , A — 1 . Moreover, by using the expression ^p) for the period T, one can simplify Eq. ([23|l to the following 
form 



Accordingly, the Floquet spectrum writes 

Afe + iuk 



^l + ^ln[l-i^(^.)]-^ 



(25) 



(26) 



For even values of A^ and for the maximal phase = vr , an explicit solution for the Floquet exponent can be 

obtained: 



A 



l+N/2 



-1 + - In 
T 



1 



a-l + 2l3^Tg (1 + e^f^^) {e^l^^ - 2ePT + ^-PT)-^ 



(27) 




FIG. 5: (Color online) Floquet spectrum A((^) as a function of the phase <^ in the case of vanishing pulse-width. Red circles 
indicate the numerical results obtained for the linear stability analysis without any approximation for A'^ = 1,000, the black 
line refers to the first-order approximate expression (|26p and the dashed blue line to the approximation (|24p . In the inset we 
magnify the region of small values of Lp. The parameters are a — 1.3, g = —1.2, and /3 = 1.0. 



C. Isolated eigenvalues 



The main consequence of the stability analysis discussed at the beginning of this section is that there exist also 
eigenvalues whose time scale is of the order a^^ ^ N^^. These exponents are those labelled with the indices N and 
iV + 1 in Eq. pT|) . The analysis of these two exponents becomes particularly relevant when they approach and cross 
the unitary circle; their analytical expression can be easily derived by investigating the case ||/LtAr_Ar+i — 1|| ~ 1 , which 
leads to a divergence of the maximal Floquet exponent, Am N . The two corresponding eigenvalues turn out to 
satisfy the equation 



^N,N+l o I 1 1 



1 - 



l±\ 1- 



4(1 - a - gE) 



2(1 - a - gE) \ V " I^^Tg 
Whenever one of such solutions become positive, this will give a leading Floquet exponent Am ^ N . 



(28) 



D. Phase diagram 



The main difference with respect to the case of finite pulse-width is that the spectral component associated with SWs 
does not reduce to a fully degenerate zero Floquet exponent. On the contrary, it actively contributes to determining 
the stability properties of the network. While the splay state is found to be unstable for finite values of f3 and 
excitatory coupling (</ > 0), a more interesting scenario is found for the inhibitory case {g < 0) . 

Let us illustrate such situations by analyzing the stability spectrum obtained for a = 1.3, /3 = 1.0 and g = —1.2. 
In Fig. [5] we show that the theoretical expression ((25)1 reproduces most of the spectrum obtained by the numerical 
diagonalization of the exact Jacobian for a network of 1000 neurons. In particular, the agreement is very good around 
the top part of the Floquet spectrum, which accounts for the stability of the network. More precisely, the least stable 
mode is that one characterized by the highest frequency (the 7r-mode), i.e. it corresponds to an up-down behavior of 
the perturbation, when passing from one to the next neuron. As this condition holds true for arbitrary values of TV, 
we are led to conclude that one cannot perfom the continuum limit straightforwardly, since this would remove those 
fluctuations that are most relevant for the stability of the network. 

Notice also that the only region where Eq. (pS]) fails to reproduce the true spectrum is close to vanishing values of 
the frequency ip. Analogously to the previous case, in this region (see the inset in Fig. [5]) one has to invoke Eq. (j^H) 
to account for the dip centered around ~ 0. 

If one has to determine the entire spectrum of a finite network, the two components arising from Eqs. (|24l26p have 
to be properly matched. After selecting equispaced modes according to the system size (the spacing being 27r/7V), 
one has to identify the value of k where the distance between the two spectral components is minimal. Below this 
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FIG. 6: (Color online) Maximal Floquet exponent Am as a function of /3 for a = 1.3 and g = —1.2. The line bordering the 
shaded region corresponds to the numerically obtained maximum in a network of = 500 neurons; red circles correspond to the 
second eigenvalue in the same network, while the solid line corresponds to the analytical expression (|27|l for A^. In particular, 
the lower border of the dark shaded region corresponds to A,r. In the inset, the numerically computed Floquet exponent (open 
red circles) is compared to the analytical expression (|28p for Ajv,jv+i. 
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FIG. 7; (Color online) Finite-frequency Floquet spectra for j3 — 0.01 (black solid line), 0.02 (red dashed line), 0.04 (green 
dot-dashed line) and 0.08 (blue dotted line) and a — 1.3, g — —1.2. 

value the component to be considered is that of LWs (see Eq. (HH)), while above such a value the spectrum is weU 
approximated by Eq. ([2^ . 

It is also instructive to investigate the dependence of the spectrum on the parameter (3. In Fig. [6] we have plotted the 
maximum Floquet exponent Xm as a function of /3. The maximum has been determined by diagonalizing numerically 
the Jacobian for a network of = 500 neurons: it corresponds to the border of the shaded region (for the sake of 
clarity, the peak has been cut out). The comparison with the results obtained from Eq. ([^^ confirms that SW modes 
account for the stability of the entire network in a quite wide range of /3 values, except for an intermediate region, 
where the maximum Floquet exponent is well reproduced by Eq. pS)) (as it can be appreciated by looking at the inset 
of Fig. El). 

Outside this region, X„i is well reproduced by the analytical expression Eq. ([27| . except for /3 < 0.125. Indeed, 
for small /3, the maximum of the spectrum occurs for \ip\ < tt, as it can be seen in Fig. [T] Moreover, the maximal 
Floquet exponent X„i remains finite even for (3 —^ 0. Finally, inside the intermediate region dominated by Xn,n+i, 
the analytical prediction ([?7|) is extremely close to the 2nd (numerically computed) exponent. Therefore, the theory 
exposed in this section predicts correctly the stability of the network, although one must carefully identify the spectral 
component that is responsible for the relevant contribution. 

There is another important conclusion that can be drawn from Fig. El by comparing two limit cases. On the one 
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FIG. 8: (color online) Phase diagram for the stability of the splay state for a — 1.3 and inhibitory coupling in the sharp spike 
limit. The solid black horizontal line divides the stable from the unstable region, while the dashed red line encircle the strongly 
unstable phase where the maximum Floquet exponent is proportional to A''. 



hand, the limit /? — > oo corresponds to (5-like pulse. In fact, one can see that the maximum Floquet exponent coincides 
with the expression derived in Ref. On the other hand, for (3 — 0, the model converges to the a ^ oo limit of the 
model considered in [2l|. In the two limits. 



6[{a - 1)T + g]{aT + g) 



(29) 



Altogether, the /3-dependence resolves the contradiction mentioned in the beginning of this section, as we can 
conclude that the pulse- width does not only affect the quantitative value of the Floquet exponent, but contributes 
also to determining the qualitative stability properties. It is therefore useful to think about the meaning of /?. It 
expresses the pulse bandwidth 1/a in terms of units. By remembering that the interspike interval is T/N, it 

is convenient to introduce the adimensional parameter 



T/N 
1/a 



(3T . 



This parameter, being the ratio between the interspike interval and the pulse- width, is susceptible of being determined 
in general contexts that go beyond the specific pulse-shape choice and model adopted in the present paper. 

Actually, the phase diagram reported in Fig. [51 confirms that r = (3T is a meaningful parameter, since stable and 
unstable phases are separated by a critical line, where r is constant independently of g. By solving = 0, we find 
that the critical value of r discriminating between stable and unstable phase is given by the following implicit equation 



2(r2 + l)e 



3r 



1 = 



(30) 



Moreover, we see that the strongly unstable regime seen in Fig. [S] arises only for a sufficiently strong inhibitory 
coupling (g < —1). 

We conclude this section by coming back to the possible reasons for the failure of the mean field approach. In the 
upper panel of Fig.[9l we show the time evolution of the field E{t) for increasing values of N and fixed a — 120. The 
oscillations around the mean value tend to decrease and this indeed suggests that it is meaningful to introduce a sort 
of average flux of spiking neurons as done in the mean field approach. One can also see that since E{t) is increasingly 
flat (for increasing N) , it is natural to expect that the stability of the splay state is very weak: The neuron-neuron 
interactions are in fact mediated by the common fleld E. It might be tempting to reduce the stability of the splay 
state to that of the single neuron exposed to a given fleld dynamics E{t), but this would amount to no more than a 
crude approximation. In other words, there is no shortcut for performing the analysis carried out in this paper. 

In the lower panel, of Fig.O the fleld evolution is reported for the same values of N and (3 = 0.14. In this case, fleld 
oscillations become faster but their size does not decrease upon increasing A''. We consider this as the major reason 
for the failure of the mean field approach. Furthermore the finite amplitude of the oscillations is also responsible for 
mantaining a flnite stability even in the limit N ^ oo . 




FIG. 9: (Color online) Time evolution of the field E{t) associated to a splay state for increasing A'' values: A'^ = 100 (black 
solid lines), 200 (blue dashed lines), 300 (red dotted lines). The upper panel refers to a fixed value of a = 120, while the lower 
panel to a constant /3 = 0.14. The other parameters are a = 1.3 and g = —1.2 



V. CONCLUSIONS AND PERSPECTIVES 



In this paper we have shown that the stability of splay states can be addressed by reducing a model of globally 
coupled differential equations to suitable event-driven maps which relate the internal configuration at two consecutive 
spike-emissions. The analytical investigation of the Jacobian in the large N limit reveals that the spectrum of 
eigenvalues is made of three components: i) long wavelengths cigcnmodes that emerge also from a mean- field approach 
a) short wavelengths; in) isolated eigenvalues, which signal the existence of strong instabilities, i.e. eigenvalues 
that are proportional to the network size. Altogether, we find that drastically different results can be found for 
different values of the ratio r between the interspike interval and the pulse-width. This leads us to conclude that the 
stability of large networks of neurons coupled via narrow pulses, crucially depends on the parameter r, thus suggesting 
that the dynamycal stability of these models demands a more refined treatment than mean-field. It will be interesting 
to investigate the role of r in more general contexts, e.g., by considering different pulse shapes (possibly adding the 
effect of delay) and different force fields (possibly including further degrees of freedom, as it naturally appears in the 
Hodgkin-Huxley model). 

Another issue that it will be worth analyzing concerns the exact stability in large but finite networks in the presence 
of finite pulse- widths. Our analysis has revealed that the first order approximations fails to reproduce even the sign 
of the maximum Floquet exponent. That is a consequence of the decay of the spectrum which, in turn, requires 

developing a third-order perturbation theory. Our numerics shows that the splay state is always stable towards 
SW perturbations in networks of LIF neurons with excitatory coupling, while a 1st and 2nd order approximation 
theories may give rise to unstable states. Is this an indication of a systematic drawback of discrete-time models, or 
an indication that qualitatively different stability properties may be found for different force fields? The analysis of 
nonlinear force fields is definitely needed for obtaining a convincing answer to this question. 
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